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O ' ABSTRACT 

<Si Aims. We present global 3D MHD simulations of disks of gas and solids, aiming at developing models that can be used to study 
various scenarios of planet formation and planet-disk interaction in turbulent accretion disks. 

Methods. We employ the PENCIL CODE, a 3D high-order finite-difference MHD code using Cartesian coordinates. We solve the 
, equations of ideal MHD with a local isothermal equation of state. Planets and stars are treated as particles evolved with an 

■ N-body scheme. Solid boulders are treated as individual superparticles that couple to the gas through a drag force that is linear 
r^*) in the local relative velocity between gas and particle. 

O'N ■ Results. We find that Cartesian grids are well-suited for accretion disk problems. The disk-in-a-box models based on Cartesian 
' grids presented here develop and sustain MHD turbulence, in good agreement with published results achieved with cylindrical 

\l codes.We investigate the dependence of the magnetorotational instability on disk scale height, finding evidence that the turbu- 
lence generated by the magnetorotational instability grows with thermal pressure. The turbulent stresses depend on the thermal 
pressure obeying a power law of 0.24 ± 0.03, compatible with the value of 0.25 found in shearing box calculations. The ratio of 

t , Maxwell to Reynolds stresses decreases with increasing temperature, dropping from 5 to 1 when the sound speed was raised by 

■ a factor 4, maintaing the same field strength. We also study the dynamics of solid boulders in the hydromagnetic turbulence, 
by making use of 10 6 La grangian particles embedded in the Eulerian grid. The effective diffusion provided by the turbulence 
prevents settling of the solids in a infinitesimally thin layer, forming instead a layer of solids of finite vertical thickness. The 
measured scale height of this diffusion-supported layer of solids implies turbulent vertical diffusion coefficients with globally 

5-j averaged Schmidt numbers of 1.0±0.2 for a model with a w 10~ 3 and 0.78±0.06 for a model with a w 10 _1 . That is, the vertical 
turbulent diffusion acting on the solids phase is comparable to the turbulent viscosity acting on the gas phase. The average 
bulk density of solids in the turbulent flow is quite low (pp=6.0 x 10 -11 kgm~ 3 ), but in the high pressure regions, significant 
overdensities are observed, where the solid-to-gas ratio reached values as great as 85, corresponding to 4 orders of magnitude 
higher than the initial interstellar value of 0.01 
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1. Introduction 

Planets have long been believed to form in disks of gas 
and dust around young stars (Kant 1755, Laplace 1796), 
interacting with their surroundings via a set of complex 
and highly nonlinear processes. In the core accretion sce- 
nario for giant planet formation (Mizuno 1980), dust co- 
agulates first into km-sized icy and rocky planetesimals 
(Safronov 1969, Goldreich & Ward 1973, Youdin & Shu 
2002) that further collide, forming progressively larger 
solid bodies that eventually give rise to cores of several 
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Earth masses. If a critical mass is attained, these cores be- 
come gas giant planets by undergoing runaway accretion 
of gas (Pollack et al. 1996). Otherwise, just a small amount 
of nebular gas is retained by the core, which ends up as 
an ice giant. 

The success of this picture in explaining the over- 
all shape of the solar system was shaken by the discov- 
ery of the extra-solar planets. In less than a decade, the 
zoo of planetary objects received exotic members such 
as close-in Hot Jupiters (Mayor & Queloz 1995), pul- 
sar planets (Wolszczan & Frail 1992), highly eccentric gi- 
ants (Marcy & Butler 1996), free-floating planets (Lucas & 
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Roche 2000), and super-Earths (Rivera et al. 2005). Thus, 
understanding the diversity of these extra-solar planets 
is a crucial task in planet formation theory. 

Planet-disk interaction seems to be one of the obvi- 
ous candidates to account for this diversity. Planets ex- 
change angular momentum with the disk, leading to ei- 
ther inward or outward migration (Ward 1981; Lin and 
Papaloizou 1986; Ward & Hourigan 1989; Masset et al. 
2006). An understanding of the physical state of accre- 
tion disks is essential to provide a detailed picture of the 
effect of migration on planetary orbits. 

Analytical theory must necessarily contain a number 
of linearizing simplifications. Therefore, numerical simu- 
lations are a major tool to provide advances in the prob- 
lem. But even then, the large computational demands of 
such calculations have put some restrictions and limi- 
tations in the models presented so far. Because of this, 
although many of the individual physical processes oc- 
curring on circumstellar environments are understood in 
some detail, state-of-the-art calculations on planet forma- 
tion still lag behind our current understanding, contain- 
ing simplifying assumptions needed to reduce the com- 
putational effort. 

For example, the evolution of temperature is usu- 
ally neglected in solving the dynamical equations, fa- 
voring an imposed temperature profile. Paardekooper & 
Mellema (2006) showed that in non-isothermal disks, the 
net torques acting on a forming planet can change sign 
due to asymmetric heating on the planet's corotation re- 
gion, potentially stopping and reversing the migration 
of the planet. 2D and 3D models of disks with radiative 
transfer were presented by D'Angelo et al. (2003) and 
Klahr & Kley (2006), showing that a high-mass planet 
may carve a cold gap in the disk while retaining a thick 
circumplanetary cloud. But no radiative global simula- 
tion with explicit ray tracing, able to consistently treat op- 
tically thin and thick regions and the transition between 
them, has been presented so far. 

Magnetic fields have been shown to play a major 
role in the structure and evolution of accretion disks. 
Observational efforts in the detection and analysis of pro- 
toplanetary disks show evidence that these disks accrete, 
with a mass accretion rate of the order w 10 -8 MQyr -1 
(e.g., Sicilia-Aguilar et al. 2004). Such a powerful ac- 
cretion cannot be explained by molecular viscosity, re- 
quiring some other mechanism to transport angular 
momentum outward. Balbus & Hawley (1991) pointed 
out the importance of the magnetorotational instabil- 
ity (Velikhov 1959, Chandrasekhar 1960, 1961) for accre- 
tion disks. In their important work, they show that this 
magnetorotational instability (MRI) is operative in suffi- 
ciently ionized Keplerian disks as long as the magnetic 
field is subthermal, generating a turbulence powerful 
enough to explain observed accretion rates in protoplan- 
etary disks. 

However, although magnetic fields are ubiquitous 
in the universe, protoplanetary disks are thought to be 
"cold" and thus not completely ionized. Cosmic rays can 



provide the required ionization for the MRI to operate, 
but they cannot penetrate all the way to the midplane of 
the disk (a standard value for the penetration depth is a 
gas column density of Z = 100 gcm~ 2 ). The result is that 
in the region where giant planets are thought to form, 
only the surface of the disk is sufficiently ionized for the 
MRI to grow. Turbulence thus likely operates in a surface 
layer, while the midplane is neutral and laminar, consti- 
tuting a so called "dead zone" (Gammie 1996, Miller & 
Stone 2000, Oishi et al. 2007). 

As a result of the mentioned difficulties of model- 
ing the coupled interaction between radiation, magnetic 
fields, dust grains, solids, neutral and ionized gas in the 
gravitational potential of a star and embedded planets, 
the numerical works in the field show a heterogeneity 
of methods, with most works tackling only some aspects 
of the problem. Particularly, numerical simulations have 
focused on local Cartesian shearing boxes (e.g., Hawley 
et al. 1995, Brandenburg et al. 1995) or global disks on 
cylindrical grids (e.g., Hawley 2001, Armitage et al. 2001, 
Nelson 2005). As the MRI is a local process, the shearing 
box has the advantage of capturing much of the physics 
of the problem while significantly reducing the computa- 
tional effort and complexity - for instance (shear-) peri- 
odic boundary conditions can be used. The global disks 
on cylindrical grids offer the advantage of having the 
grid and flow geometry coinciding, but in this case, spe- 
cial care must be taken for the boundary conditions, as 
reflective boundaries make waves bounce through the 
computational domain in an unphysical manner and out- 
flow boundary conditions may lead to too much mass 
loss (Fromang & Nelson, 2006). Fromang & Nelson (2006) 
have also presented the first simulation of the MRI in 
global disks with vertical density stratification. A com- 
parison between their models and a stratified version of 
ours will be addressed in future work. 

In a series of articles we aim at constructing global ra- 
diative magnetohydrodynamical simulations. In this first 
paper, we present the features and capabilities of the 
numerical scheme used by constructing cylindrical disk 
models of gas and solids with MHD turbulence. These 
models will be developed in future work to allow for 
stratification, radiation and a global self-consistent treat- 
ment of dead zones. 

The simulations presented here are embedded in 
Cartesian boxes. Although it can be regarded as unpracti- 
cal for simulating a flow with cylindrical symmetry, such 
a grid also presents some advantages. First, cylindrical 
grids are a strong limitation for flows that deviate from 
cylindrical symmetry, e.g. circumbinary disks or 3D jet 
simulations, mainly because it is impossible to have a 
flow across the center of the grid, and at r=0 reflection 
must occur. Second, this approach has proved useful in 
view of computational simplicity and parallelization ef- 
ficiency (e.g. Dobler et al. 2006). In particular, by hav- 
ing cells of constant aspect ratio, Cartesian grids have 
much reduced numerical dissipation when compared to 
grids with complex geometry (van Noort et al. 2002). 
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Therefore, while cylindrical and spherical grids can ex- 
plicitly conserve angular momentum w.r.t the origin of 
the coordinate system, it is of no benefit for systems that 
do not have the center of mass at the origin. Third, as 
photons travel in straight lines (in the absence of gen- 
eral relativistic effects), a radiative transfer scheme with 
ray tracing is simpler to implement in a Cartesian grid in 
spite of the cylindrical symmetry of the hydrodynamical 
flow (Freytag et al. 2002). 

As a numerical solver we employ the PENCIL CODlfl, 
a high (6th) order finite difference code. Such a numeri- 
cal tool is highly different from most other astrophysical 
codes in use in the literature (see de Val-Borro et al. 2006 
and references therein), thus also providing an indepen- 
dent check of the results so far obtained in the field. 

This paper is structured as follows: we discuss the 
model in Sect. 2, proceeding to test cases in Sect. 3. In Sect. 
4 we discuss the several MHD simulations performed. In 
Sect. 5 the models with solids are presented, finally lead- 
ing to the conclusions in Sect. 6. 

2. The model 

2.1. Gas dynamics 

The equations solved are those of ideal MHD in an in- 
ertial reference frame with a central gravity source. The 
equation governing the evolution of density is the conti- 
nuity equation 



§f = -pV-u + f D (p), 



(1) 



where p and u are the density and velocity of the gas. The 
operator D/Dt = d/dt + u ■ V represents the advective 
derivative. 

The equation of motion is the sum of all forces acting 
on a parcel of gas. It reads 

/ x B 



Du 1„ 
— = --Vp 

Dt p v 



V<3> 



fv{u,p), 



(2) 



where p is pressure, <J> the gravitational potential, B is the 
magnetic field, / = fi^ 1 V x B is the volume current den- 
sity, as defined by Ampere's Law, and }Iq is the magnetic 
permeability of vacuum. 

The evolution of the magnetic field is governed by the 
induction equation. The PENCIL CODE, however, works 
not with the magnetic field itself, but with the magnetic 
potential A, where B = V x A. This automatically guar- 
antees the solenoidality of the magnetic field, as the con- 
dition V-B = V- (VxA)=0is always satisfied. The 
induction equation formulated for the magnetic potential 
reads 

3A 

— =uxB + f v (A). (3) 

The equation of state, relating pressure and density, 
closes the system of equations. We use the ideal gas law 



(4) 



with a locally isothermal approximation, where the 
sound speed c s is a time-independent function of the 
cylindrical distance s to the z-axis. We write cylindrical 
coordinates as (s,(p,z) and spherical coordinates as (r,(p,6), 
where 8 is the polar angle and (p the azimuthal angle. The 
z direction is perpendicular to the midplane of the disk. 

The gravitational potential <1> has contributions from 
the star and the N — 1 embedded planets, 



GMj 



n 
V 



(5) 



where G is the gravitational constant, M, is the mass of 
particle i and 72.; = \r — r Pi | is the distance of a gas parcel 
relative to particle i. The quantity b\ is the distance over 
which the gravity field of the particle i is softened to pre- 
vent singularities. 

The functions f D (p), f v (u,p), and f,j(A) are explicit 
mass diffusion, viscosity and resistivity terms, needed to 
stabilize the numerical scheme. They are composed of 
two terms, where the first one is a conservative sixth- 
order dissipation. This term is described in detail in 
Haugen & Brandenburg (2004) as well as in Johansen 
& Klahr (2005) for the case of isotropic dissipation. A 
generalization for the anisotropic case, required for non- 
cubic cells, is shown in Appendix [Aj The second term is 
a localized shock-capturing dissipation, activated when 
large negative divergences, typical of shocks, are formed 
(Haugen et al. 2004). This is described in Appendix [B] 

2.2. Planet orbital evolution 

The star and the planets are treated as an N-body ensem- 
ble, evolving due to their mutual gravitational interac- 
tion. The equation of motion for particle i is 



dv Pi 
~~dt 



£ GM i<* 



(6) 



where 1Z,j = |r p; — r Pj \ is the distance between particles i 
and and 7c;j is the unit vector pointing from particle / 
to particle i. The first term of the R.H.S. is the combined 
gravity of the gas onto the particle i 



{Tl\ + b\) 



3/2 



dV, 



(7) 



See http:/ /www.nordita.dk/ software /pencil-code 



where the integration is carried out over the whole disk. 
As we are not interested in the disk's self-gravity, but 
rather on its gravitational effect on one specific point (or 
a few points in case of multiple planets), calculating the 
integral above is simpler and faster than using a Poisson 
solver to find the gravitational potential of the disk ev- 
erywhere on the grid. 

The smoothing distance bj is taken to be as small as 
possible. It is usually a fraction of the Hill radius. For 
reasons described in sect. 2.7, the stellar potential can be 
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treated as unsoftened (£v=0). We note that in this formu- 
lation there is no distinction between a planet and a star 
except for the mass. The star evolves dynamically due to 
the gravity of the planets, wobbling around the center of 
mass of the system, which is set to the center of the grid. 
As the disk is not massive compared to the star, we ex- 
clude the disk torques from influencing the star. For runs 
without planets a constant gravity profile with a star at 
the center of the grid is used instead of solving the equa- 
tions of the N-body code. 

2.3. Dynamics of solids 

To model the early stages of planet formation where 
solids grow from cm and m sizes to kilometer-sized plan- 
etesimals we consider the dynamics of meter-sized solid 
boulders, also treated as individual Lagrangian particles. 
Each of the particles has its own position and velocity, 
independent of the grid, integrated by the same particle 
module of the PENCIL CODE that is used for the plan- 
ets. The difference is that as the planets interact with the 
disk and with themselves by gravity, the particles inter- 
act with the disk only via a drag force that is proportional 
to the velocity of the particle with respect to the local gas 
velocity. 

While in our cylindrical models there is no vertical 
gravity on the gas, the particles do feel this component 
without which no settling towards the disk midplane 
would occur. The evolution equation for solid particle i 
is therefore 



where iy is the friction time and u is the gas velocity 
at the position of a particle. We assume that the fric- 
tion time is independent of velocity differences between 
gas and particles. We choose it to be T=l/Qo, which 
for the typical densities and temperatures in the disk 
(Table 1), corresponds to particle radii between 0.4 and 
2.5 meters, depending inversely on the orbital distance. 
The assumption of linearity of the drag law holds as 
long as the velocity difference between gas and solids 
is much smaller than the sound speed (Weidenschilling 
1977, Paardekooper 2007, Johansen et al. 2007). The con- 
dition is met since the turbulence generated by the MRI 
is subsonic. 

The gas velocity u at the position of the particle is 
interpolated from the nearest 27 grid points, using a 
Triangular Shaped Cloud scheme, as described in Youdin 
& Johansen (2007). 

2.4. The code 

The PENCIL CODE is a non-conservative Cartesian finite- 
difference MHD code that uses sixth order centered spa- 
tial derivatives and a third order Runge-Kutta time- 
stepping scheme, being primarily designed for compress- 
ible turbulent hydromagnetic flows. 



Table 1. Conversion between code and physical units 



Quantity 


Physical Unit 


Length 


5.2 AU (= 7.8 x 10 n m) 


Density 


2.0 x 1CT 8 kgm~ 3 


Velocity 


L.jL x 1U ms 


Energy 


1.60 x 10 36 J 


Pressure, Stress 


3.41 Pa 


Time 


1.89 yr (= 6.0 x 10 7 s) 


Magnetic Field 


2.07 x 10~ 3 T 


Viscosity 


1.02 x 10 16 m 2 s -1 


Mass 


4.73 x lO- 3 M e 


Mass accretion rate 


2.51 x 10~ 3 M e yr" 1 


Domain Size (L S ,L Z ) 


2-13 AU, ±1.3 AU 


Resolution (A x) 


0.08 AU 



The PENCIL CODE was recently applied to a 2D global 
laminar disk calculation, in which the results agreed with 
those of polar-grid based codes (de Val-Borro et al. 2006). 
We extend this calculation now to three dimensions with 
magnetic fields, fully exploiting the capabilities of the 
PENCIL CODE for handling the problem of numerical hy- 
dromagnetic turbulence. 

2.5. Units 

We adopt dimensionless units such that 

GM = po = no = 1 

The quantity GM has dimension of length 3 time -2 , so it 
sets a constraint on [x][f]. The unit of time follows from 
this as being the inverse of the Keplerian angular fre- 
quency at s = Sg = 1 

w = /f ="„-', » 

which gives an orbital period P = 2tt at So in absence of 
a global pressure gradient. 
The unit of velocity 

[u] = [x]/[t] = n s 

is therefore the local Keplerian speed at So- The sound 
speed is set accordingly, through the Mach number (see 
equation fTT|l . Density is measured relative to the initial 
density of the box [p] = Pq. 

The unit of magnetic field follows from the Alfven 
speed, 

[B] = OoSoVWo- 

It follows from this that the unit of magnetic vector po- 
tential is 

[A] = [B][x] = QqSoVWo- 
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As the simulation is dimensionless, it scales with the 
choice of physical units. By assuming that Sq is the semi- 
major axis of Jupiter, flj = 5.2AU, and considering the 
typical density of the minimum mass solar nebula at that 
location, po w 2 x 10~ 8 kgm -3 , the physical units corre- 
sponding to the employed code units are listed in Table 1. 



2.6. Initial conditions 

We use a Cartesian box with a spatial range x, y £ 
[-2.6,2.6], and z £ [-0.26,0.26] (see Table 1 for a con- 
version of the units used to physical units). The num- 
ber of cells is usually N x =Ny=320, N z =32. This ensures 
that A x= A y=A z, i.e, all cells are cubes of the same size. 
However, for some models we double the resolution in 
the vertical direction in order to resolve faster growing 
wavelengths of the MRI. Doubling the resolution in x and 
y would keep the cells cubic, but the already expensive 
computational costs would become unpractical without 
yielding any other major advantage. We therefore keep it 
at 320 x 320 and introduce anisotropic hyperdif f usivity to 
treat the non-cubic cells (see Appendix |Al. 

As stated before, we use the ideal gas law approxima- 
tion to evaluate the pressure. The sound speed is set as a 
power law 



Border Profile 



-q T /2 



c So s ->t'~. (10) 

We usually set cj T =l, so that the Keplerian flow has a con- 
stant Mach number Ai 



M 



Cl K s 



const., 



(11) 



where is the "cylindrical" Keplerian angular velocity 
profile 

n2 _ GM* 

L1 K-—- 

The Mach number is seen to be the inverse of the aspect 
ratio h = H/r, where H = c s /Ci is the pressure scale 
height. We checked the evolution and saturated state of 
the turbulence for Mach numbers of 5, 10 and 20. 

We also perform simulations with radially-varying 
A4 where the sound speed follows a steeper power law, 
with q T =2. 

Accretion disks exhibit a radial density gradient, but 
this gradient arises due to accretion itself (Shakura & 
Sunyaev 1973). Therefore, we initialize the midplane den- 
sity at a constant value po, in order to understand the role 
of the stresses in generating the density gradient. 

We start our models in strict equilibrium between 
gravity, global (thermal and magnetic) pressure gradients 
and centrifugal forces, 

d ( ]?_ 



Q = Cli 



sp ds 1 



(12) 



2. 7. Boundary conditions 



In this work, we compute models with and without an 
inner boundary to quantify the advantages /drawbacks 



1.2 
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0.4 
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Fig. 1. Border profile of the simulations. The step func- 
tion (Eq. [131 is applied to the derivatives of the dynam- 
ical variables. The gas is free to evolve between s=0.5 
and s=2.3, is slowed down between [0.4,0.5] and [2.3,2.5] 
(light-shaded areas) and is effectively frozen at s < 0.4 
and s > 2.5 (dark-shaded areas). 

of such a feature in a Cartesian grid. As for the external 
boundary, the box limits at x = ±2.6 and y = ±2.6 do 
not correspond to the physical boundaries of the prob- 
lem. Indeed, the dynamically evolving disk encompasses 
a cylinder inside of s ex t = 2.5. The frozen regions outside 
of this cylinder play the role of ghost rings in cylindrical 
codes. 

After evolving the dynamical equations, we set the 
time derivatives of all variables to zero in the region out- 
side s ex t = 2.5. As the variables cannot evolve outward 
of s ex t , being effectively frozen in this region, the "real" 
boundary conditions of the box (e.g., open, reflecting) do 
not matter if this freezing boundary condition is used. 

To avoid numerical instabilities due to this abrupt 
jump from frozen to evolving regions, we apply a buffer 
zone to the derivatives of the variables, that smoothly 
drives the variable X to a desired value Xo in a timescale 
T, such that cp 



dX 
Ht 



X-X 



S(s). 



(13) 



where S(s) is the uniquely defined fifth-order step func- 
tion that is 1 at the domain boundary and at the interior 
boundary of the buffer zone while maintaining continu- 
ous second order derivatives (Dobler, private communi- 
cation). Its shape is visualized in Fig. 1. We usually take 
the driving term Xo to be the initial condition of the vari- 
able, and the driving time T being the Keplerian period 
In I Qk/ the dynamical timescale of the disk. The effect of 
this border profile is to smooth the transition between the 
evolving disk and the frozen regions of the grid, thus pre- 
venting large gradients and discontinuities that would 
otherwise arise. As the gas flow is symmetric in the ver- 
tical direction, the vertical boundary condition is set to 
periodic for the purpose of simplicity. 
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For the runs without an inner boundary, we smooth 
the quantities containing singularities by replacing 

s -" => (s 2 + b 2 )-' l/2 . (14) 

In practice, it is applied to the angular frequency Cl, the 
gravitational potential <3> and the sound speed c s . We usu- 
ally take b = 0.1, so the smoothed gravitational potential 
deviates from the Newtonian by less than 5% at s^t = 0.4. 
The physical domain thus runs from s^t to s e xt • 

For runs with an inner boundary, we apply inside 
Sjnt the same freezing as used outside s ex t . The N-body 
particle code does not participate in the freezing, so al- 
though the star lies in a region of frozen gas, it is allowed 
to move. 

As the gas is frozen in the inner and outer parts, the 
information about the flow in this region is not of interest. 
Therefore, we exclude these regions from the time-step 
calculation. As their time derivatives are set to zero at the 
end of the time-step, they cannot violate causality. 

In principle, we could set as close to zero as pos- 
sible (by not using smoothing but retaining an inner 
boundary), in order to study the processes that hap- 
pen in the immediate vicinity of the star, like winds, 
the magnetic cavity and surface accretion (von Rekowski 
& Piskunov 2006). However, due to the increasing 
Keplerian velocity in the advection and the decreasing 
resolution of the orbits, non-axisymmetric wave modes 
(particularly the m=A mode) build up in the inner disk as 
we try to push — > 0. The density fluctuations result- 
ing from the excitation of these modes lead to numerical 
instabilities. 

Finally, the magnetic potential follows the same 
boundaries as described above. This would be a prob- 
lem if we solved for the actual magnetic field, as sinks 
or sources of magnetic flux imply the presence of open 
magnetic loops (monopoles). By solving for the magnetic 
vector potential we do not face such problems. 

The solid particles obey different boundary condi- 
tions, explained in Sect. 5. 

3. Influence of free parameters 

In order to clarify the influence of the numerical scheme 
and the approximations made, a series of non-magnetic 
2D models were computed, with and without planets. 
The grid being Cartesian, all our simulations span the 
whole azimuthal domain. We usually evolve the simu- 
lations up to 100 orbits at sq, which corresponds to «25 
orbits at the outer edge of the disk and w400 orbits at its 
inner edge. 

3.1. Viscosity 

Explicit hyperviscosity and hyper diffusion induce dissi- 
pation primarily near the grid scale, replacing the usual 
2 nd order Laplacian terms. A visual picture of the differ- 
ence between using the two types of viscosity is seen in 



Fig. 2. The first model was computed with a Laplacian 
viscosity V\ = 10~ 3 . The radial inflow is significant, and 
as the outer frozen region behaves like an infinite reser- 
voir of matter, the total mass inside the disk keeps on 
rising as matter flows in from this reservoir. The radial 
density profile soon starts to deviate from the flat initial 
condition . Shown in the figure is the density profile after 
100 orbits at sq with v\ = 10 -3 . When using hyperviscos- 
ity of similar strength at the grid scale, i.e., 1/3 = 10 -10 , 
the overall flow shows no significant deviations from the 
initial conditions. 

The simulation shown using Laplacian viscosity was 
computed without an inner boundary, using a softened 
stellar potential with b=Q.l. Using the damping and freez- 
ing profile described in Sect. 2.7, the density is not al- 
lowed to deviate much from the initial condition at the 
boundaries. In the physically evolving part of the disk, 
however, a density profile of exponent s~ 0A evolves. 

3.2. Mass diffusion 

To evaluate the influence of mass diffusion, we simulate a 
laminar disk with constant Laplacian viscosity V\ = 10~ 5 
in the presence of a gap-opening Jupiter-mass planet. 
We performed runs of resolution 320 x 320 with hyper- 
diffusion coefficients ranging from D 3 = 5 x 10~ n to 
10~ 14 . After a hundred orbits, time enough for the planet 
to open a deep gap, the density profiles are plotted in 
Fig. 2a, where a run with resolution 640 x 640 without 
explicit diffusion is shown for comparison. 

It is seen that D3 =5xl0~ n constitutes too much 
diffusion, as the gap is significantly altered. As less dif- 
fusion is used, the shape of the gap monotonically ap- 
proaches the one recovered in the higher resolution run. 
The walls of the gap are fairly well reproduced for lower 
diffusion, but its bottom is always shallower even for the 
lowest coefficient used (D3 = 10~ 14 ). 

Judging from the gap alone, one could in principle use 
no diffusion at all, but the inner disk suffers depletion for 
low diffusion regimes, due to the non-conservative na- 
ture of the numerical scheme. Even the higher resolution 
run seems to have lost mass due to the lack of explicit 
diffusion. 

We adopt a hyper diffusion coefficient of D3 = 5 x 
10 -12 as the best compromise between the need for pre- 
serving material in the inner disk and for reproducing the 
overall shape of the gap. Requiring Schmidt and mag- 
netic Prandtl numbers of 1 at the grid scale, we set hy- 
perviscosity and hyperresistivity to the same value. 

3.3. Shocks 

Shock viscosity and shock diffusion are needed for 
two reasons: (a.) to stabilize the flow near the shock- 
generating particles in runs with planets and (b.) to treat 
eventual supersonic motion in the turbulence (arising 
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Fig. 2. The radial density profile after 100 orbits for 
Laplacian viscosity v\ = 10~ 3 is compared to the pro- 
file obtained by using sixth-order hyperviscosity (V3 = 
10 -10 ) of same strength in the small scales. In the hy- 
perviscous case, the global flow is unaffected. Note that 
the frozen regions behave like infinite reservoirs of mat- 
ter. The power law describing the resulting density pro- 
file for normal viscosity is very close to the s~ - 5 , as ex- 
pected for constant viscosity (see Pringle 1981 and refer- 
ences therein) 



when the disk is exposed to a strong net vertical field), 
in which case shock resistivity is also included. 

In Fig. 2b we show gap-opening runs with fixed hy- 
perdiffusion coefficient D3 but varying the shock diffu- 
sion coefficient. From the continuity equation, one can 
tell that the effect of shock diffusion is to slow down the 
time evolution of density by smearing out any large di- 
vergences. Indeed, one sees that after a hundred orbits, 
a shock diffusion coefficient of 1 fails to reproduce the 
shape of the gap as compared to the higher resolution 
run without shock diffusion, while 10 _1 shows less accu- 
mulation than expected in the Lindblad resonances, also 
seen as compared to the higher resolution run. We there- 
fore use shock diffusion of 10~ 2 for the turbulent runs. 
The flow around a high-mass planet, however, could only 
be stabilized with a shock viscosity of 1 . 

Shock resistivity is also used when the run involves 
the magnetic potential. The value used was not tuned in 
2D runs like shock diffusion and shock viscosity. Instead, 
in the presence of turbulence, we simply checked what 
was the lowest shock resistivity coefficient that did not 
lead to numerical instabilities for model A (see Table 2), 
finding that it is of the order of unity, like the shock vis- 
cosity. 

3.4. Non-turbulent runs 

To verify the numerical stability of the model in the ab- 
sence of physical turbulence, we perform tests for cases 
where the turbulence is not supposed to be present. In 
these runs, we monitor the evolution of the mass inflow 
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Fig. 3. Upper Panel. The gap carved by a 1 Mj planet in 
a 2D disk reveals the influence of explicit diffusion in 
the calculations. The inner disk loses mass depending on 
the amount of diffusion. The value of D3 = 5 x 10 -12 
seems to ensure mass conservation in the inner disk yet 
not distorting the shape of the gap. The solid line repre- 
sents a 640 x 640 run without diffusion, for comparison. 
Resolution is 320 x 320 otherwise. 

Lower Panel. Same but with hyperdiffusion set to D3 = 
10" 11 and varying the shock diffusion coefficient from 
10~ 2 to 1. 

rate M defined as the ID radially dependent surface in- 
tegral over a surface A which is a cylinder at a radial dis- 
tance s from the origin 



M(s) 



'A 
2ns 



pun dA 

l z /2 



p(s,z)u s {s,z) dz. 



(15) 
(16) 



-Lz/2 



We see that in a 2D laminar model, the mass inflow rate 
is constant through the radial domain once a steady flow 
is achieved (which simply states that mass is conserved). 
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We thus define the mass accretion rate as the mass inflow 
rate across the inner boundary, 

th = M| s=04 , 

meaning that after crossing this boundary, the matter 
is considered lost (accreted). 

We also measure the kinetic alpha parameter of tur- 
bulent viscosity, defined as 

2R S <P 

«R — o 2' 

3 pcj 



Kinetic alpha 



where — p5u s 5u,p is the Reynolds stress; and its mag- 
netic counterpart 



2M S * 



3 pc 2 s 



where = fi^SBgSBip is the Maxwell stress. 

In Fig. |4]we show a 2D run where the velocities were 
perturbed with noise of u rms = 10 -2 , but the induc- 
tion equation is not solved. The turbulence dies out so 
fast that even before ten orbits at Sq the flow is already 
smooth, showing that the numerical scheme does not 
spuriously generate or sustain turbulence. 

A 3D cylindrical run in which we add a vertical net 
field of strength Bq = 10 -3 (dimensionless), correspond- 
ing to plasma ,6=5000 at s , 12500 at s^and 2000 at 
Sext (see Eq. I1T91 ), but where the initial flow is not per- 
turbed by noise, does not develop turbulence. We also 
tested if the MRI would develop in a disk seeded only 
with noise in both the velocity and magnetic potential 
(frms = -Arms = 10 -4 ). There is a short growth in mag- 
netic energy, presumably due to reconnection of the field 
lines, but without a structured field to maintain the tur- 
bulence the Reynolds and Maxwell stresses quickly level 
down to zero as viscosity and resistivity smooth the im- 
posed noise. 

4. Cylindrical disk runs 

For the main simulations in this paper we consider flat 
vertical profiles for the gravity field. Such an approxima- 
tion is called a "cylindrical" disk and has been often used 
in order to study the MRI (e.g., Armitage 1998, Hawley 



2001). The vertical gravity 



-Q 2 z is switched off 



so that, physically, the star is no longer a point mass at 
r=0, but a rod at s=0 extending through the length of the 
z-axis. In such a setup, the pressure scale height H has 
no hydrostatic meaning, being only a way to write the 
temperature profile of the disk. We performed simula- 
tions with non-zero net flux magnetic fields B = Bqz and 
B = £>o</>. Models with a radially varying vertical field 
proportional to Q(s) were also computed. 

For the turbulence to develop, the unstable wave- 
lengths of the MRI must be resolved. The characteris- 
tic vertical wavelength Agn of the hydromagnetic turbu- 
lence is given by (Balbus & Hawley 1991, 1998) 



A 



BH 



(17) 
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Fig. 4. Without solving the dynamical equations for the 
magnetic field, even initially vigorous random motions 



Of Mr 



0.2c Sq die out rather quickly. The plots show the 



time evolution of the globally averaged Reynolds stress 
and mass accretion rate. Time is quoted in orbits at Sq. 



where v A is the Alfven speed 



(18) 



The turbulence will be present as long as the critical 
wavelength A c is resolved. This wavelength is A c = 
^BH / v3/ whilst the most unstable wavelength of the 
MRI is 4A BH /\/l5 (Balbus & Hawley, 1991). The plasma 
B parameter - the ratio of thermal to magnetic pressure 
- can be expressed in terms of the sound and Alfven 
speed, by writing the magnetic pressure Pm = B 2 / (2«o 
in terms of the Alfven speed Pm = B 2 /(2mq) 



K/2- 



2cf_ 
v 2 ' 



(19) 



The constant Bo is usually set to 10 3 , but runs varying 
the field from 10" 4 to 10 _1 were also studied. Although 
the runs reported here are in the local isothermal approx- 
imation, we also varied the initial sound speed at Sq, c SQ , 
in order to check how the resulting turbulent viscosity 
depends on the global gas pressure. 

The parameters of the cylindrical models presented 
here are specified in Table 2. 

4.1. Constant vertical field - model A 

The evolution of the turbulence in a fiducial run with a 
net vertical flux of strength Bo — 10~ 3 and temperature 
profile corresponding to c Sg = 0.05 is shown in Fig. [5] 
The absolute value of the Maxwell stress at saturation is 
always larger than the Reynolds stress, but the latter fluc- 
tuates more strongly. 

The minimum ratio of stresses is — M S< VR S< ^=3 (at 
tw75 orbits), but it reaches as much as 100 (at t«58 or- 
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Table 2. Cylindrical turbulent disk models. 
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Parameter Results 



Run 


s int 


Bo 






fa 


N z 


Xsh 




—M s<p 


OCR 




Brms 


ft 


St 




(xlO 3 ) 












(xlO 5 ) 


(XlO 5 ) 


(xlO 3 ) 


(xlO 3 ) 


(XlO 3 ) 




(XlO 3 ) 


Uniform Field B z 


A 


0.4 


1 


0.05 


1 


5000 


32 


i 


10 6 0.24±o.04 


1.5±0.3 


0.9±o.2 


6±i 


17±9 


13±3 


7±i 


B 


0.4 


1 


0.10 


1 


20000 


32 


i 


1.0±0.2 


2.5±o.3 


0.7±o.i 


1.8±0.2 


16±5 


65±7 




C 


0.4 


1 


0.20 


1 


80000 


32 


i 


4±i 


5.3±o.8 


0.9±o.2 


1.3±0.2 


26±7 


81±7 




A2 


0.0 


1 


0.05 


1 


5000 


32 


i 


0.20±o.o4 


1.3±o.i 


0.7±o.i 


4.7±o.3 


17±5 


13±i 




B2 


0.0 


1 


0.10 


1 


20000 


32 


i 


0.9±o.2 


2.6±o.3 


0.8±o.i 


2.1±0.2 


20±n 


43±i3 




C2 


0.0 


1 


0.20 


1 


80000 


32 


i 


5±3 


5±i 


1.2±0.8 


1.2±0.3 


22±9 


U6±i9 


















Radially Varying Field j 


% 












D 


0.4 


20 


0.10 


2 


12 


64 


2 


10 6 22± 2 


78±7 


25±i 


87± 3 


71 ±9 


4±i 


140±io 


E 


0.4 


20 


0.20 


2 


50 


64 


2 


35±7 


87±is 


13±3 


30±5 


60±23 


ll±i 




Dw 


0.4 


5 


0.10 


2 


750 


64 


2 


4.1±0.4 


13±2 


5.4±o.e 


17±2 


28±io 


13±2 




Ew 


0.4 


5 


0.20 


2 


3000 


64 


2 


13±4 


27±2 


4±i 


8.7±o.z 


36±n 


33 ±3 




Uniform Field By 


F 


0.0 


30 


0.05 


1 


5.5 


32 


2 


0.7±o.i 


2.9±o.4 


2.5±o.4 


ll±i 


12±4 


24± 2 




G 


0.0 


30 


0.20 


1 


90 


32 


2 


5±2 


18±6 


0.9±o.4 


3.3±i.4 


28±n 


107±22 





bits). After 75 orbits, the average ratio of Maxwell to 
Reynolds stress is around 5. 

In agreement with previous shearing box simulations 
(Brandenburg et al. 1995, Hawley et al. 1995, Johansen & 
Klahr 2005), global disks (Hawley 2001, Nelson 2005) and 
analytical calculations (Balbus & Hawley 1991), a large 
scale toroidal field is seen to form, which dominates the 
magnetic energy, being 2 orders of magnitude stronger 
than the radial and vertical fields. Indeed, one can barely 
distinguish between the energy stored in the azimuthal 
field and the total magnetic energy. The kinetic energy is 
more evenly distributed, but it is not isotropic. The ra- 
dial component accounts for 45% of the total energy, be- 
ing w 1 .5 bigger than the vertical and twice as big as the 
azimuthal. The radial structure of the alpha parameter is 
seen in Fig. [6] The outer disk is more turbulent due to the 
smaller values of plasma S when compared to the inner 
disk. 

Following the time evolution of the turbulence in the 
midplane, it is seen that different regions of the disk reach 
saturation at different times. The turbulence starts from 
the outer disk, propagating inwards. It is expected since, 
as n decreases with radius, a uniform field implies a 
Balbus-Hawley wavelength that increases with distance 
from the star. As longer wavelengths - comparable to the 
length of the box - are easily resolved, the outer disk goes 
turbulent first. It is seen that inside Sg=l the disk did not 
go turbulent. In this model, the Mach number is constant, 
.M=20, with a constant field, B=10~ 3 through the whole 
domain, corresponding to plasma 6 running from 2000 
at s ex t and 12500 at . The magnetic field determines 
the value of the critical wavelength, which ranges from 
A c =0.002 at s int to A c =0.025 at s ext • As we have 32 grid 



points in the vertical direction, the smallest wavelength 
resolved with significant accuracy (8 points) by our high- 
order finite-difference method is L z /4«0.12. From the 
dispersion relation of the MRI (Hawley & Balbus 1991), 
this unstable wavelength has a growth rate of 0.1Q, 
much lower than the fastest growing wavelength with 
co = (3/4)n. 

We also computed models where the whole disk goes 
turbulent (models D & E), but we will use these weak 
field disk models (models A to C2, see Table 2) in the next 
subsection for studying the behavior of the turbulence 
with thermal pressure. Although dissipative and slowly 
growing, the weak field used in these disks has one ma- 
jor advantage: the turbulence grows slowly and has less 
spatial variability. Therefore, the damping at the frozen 
boundaries is more gentle than in other, rapidly grow- 
ing, violently fluctuating, disks. These issues discussed 
above reflect the compromise between keeping the disk 
cold while still retaining the field subthermal and with 
sufficient resolution to resolve the rapidly growing un- 
stable wavelengths. 

The runs with a vertical net field of Bo = 10~ 4 or 
lower did not develop turbulence as the field is weaker 
than needed in the presence of the chosen dissipation pa- 
rameters (Hawley & Balbus 1991). 

4.2. Dependence on sound speed - models B and C 

We investigate the dependence of the saturated state on 
the imposed sound speed profile. We test three different 
sound speeds of, c, =0.05, 0.10 and 0.20, corresponding 
to runs A, B and C. 
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Stresses and 




t/(2%a'o) 

Fig. 5. Time evolution of the turbulence for model A (constant net flux vertical field B = 10~ 3 z and constant Mach 
number). The top panel shows the evolution of the si^-component of the Maxwell and Reynolds stresses, while the 
middle and bottom panels show the evolution of magnetic and kinetic energy, respectively. The units are given in 
Table 1. The solid lines in the two bottom panels show the total energy. The toroidal field dominates the magnetic 
energy to the point that the energy in the azimuthal component can barely be distinguished from the total energy. 
The kinetic energy is more evenly distributed among the three dimensions, but the turbulence is not isotropic. The 
radial component shows 1.5 times more energy than the vertical and 2 times more than the azimuthal. Time is quoted 
in orbits at sq. 



Alpha - Azimuthal and Vertical Average 




Fig. 6. Radial structure of the total alpha parameter a# + 
&m for model A . The turbulence starts from the outer 
disk, where plasma /3 is smaller. The different curves cor- 
respond to snapshots at 70 (dotted), 80 (dot-dashed), 90 
(dot-dot-dot-dashed), and 100 orbits (solid), after satu- 
ration is reached. The variability is not monotonic, but 
highly fluctuating. 



As the runs are locally isothermal, we are mainly 
investigating how the MRI responds to different radial 
pressure gradients. The cold model (A) shows a weaker 
turbulence at saturation than the hotter ones, as shown 
by the strength of the stresses (Fig. [7h,b) and mag- 
netic/kinetic energies (Fig. Et,d). The Maxwell stress is 
three times bigger for model C than for A, the colder ver- 
sion. Such behavior was reported by Sano et al. (2004), 
indicating a power law of exponent 0.25 for the growth 
of the Maxwell stress with gas pressure. Our global disk 
calculations agree well with this value. 

The dimensionless magnetic parameter, which is a 
measure of turbulent viscosity, decreases drastically with 
sound speed (Fig.[8]|. As seen before, the stresses actually 
increase with increasing temperature, so this decrease of 
ocm is due to the stresses increasing less rapidly than the 
temperature. Even though alpha decreases, the effective 
viscosity v t =ctc s H increases. As seen in Fig. and Fig.[T0l 
a centrally concentrated density profile has developed 
from the initially flat configuration, a signature of mass 
accretion due to turbulent angular momentum transport, 
as also confirmed by the measured stresses. The result- 
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ing density profile in model C is smoother overall, but 
the overdensities seem to be similar in average. 

It is clear that alpha per se is not a good measure of vis- 
cosity. Since V\ = occ^ Q -1 , with no reference to the Alfven 
speed, the resulting alpha value of turbulent disks where 
the origin of the turbulence is magnetic may change with 
sound speed. As most of the analyses of turbulent thin 
accretion disks have focused on locally isothermal simu- 
lations using c s w 0.05, such dependence of a on sound 
speed did not receive proper attention. Although proto- 
planetary disks are thin, this rise in angular momentum 
transport with temperature suggests that the effects of ra- 
diation will be important for high temperature regions 
around forming planets as well as regions where the tur- 
bulence leads to significant Joule and/ or viscous heating. 

4.3. Excluding the inner boundary 

The correct treatment of boundaries is a major issue for 
numerical simulations. In global simulations of disks, 
outflow or frozen boundaries are usually used, both 
of them being more realistic than reflective, but also 
presenting disadvantages. Strictly speaking, a "perfect" 
boundary might well not exist. The best solution would 
be, of course, not having to use a boundary at all. 

Cylindrical grids only make sense if the center of 
mass is at the origin of the coordinate system as only 
then angular momentum transport can be written in a 
conservative form. If the center of mass is not at the ori- 
gin or there is more than one massive object the cylindri- 
cal coordinate system artificially introduces a reflective 
boundary at the center. Cartesian grids, as stated before, 
are not hindered by this, and we can therefore study how 
the presence or absence of an inner boundary affects the 
results. 

We compute a version of model A where the compu- 
tational domain extends all the way to Sj n t =0. Without an 
inner boundary, the gravitational potential, the angular 
frequency and the sound speed have to be smoothed ac- 
cording to Eq. ((Till to prevent singularities. When taking 
global averages, we exclude the smoothed region. 

The evolution of the turbulence and the globally av- 
eraged stresses at saturation are almost identical to those 
seen at model A. The only noticeable difference is that 
the highly fluctuating magnetic field observed in the in- 
ner boundary of model A (with rms amplitude « 2 times 
seen in the rest of the disk) does not occur in this model. 

In model A, the magnetic potential at the inner 
boundary remains frozen at the initial condition A=0. As 
the MRI builds up the magnetic potential in the freely 
evolving disk, a sharp radial gradient appears at the 
boundaries. This sharp radial gradient in the azimuthal 
and vertical components of the magnetic potential trans- 
lates into high values of the magnetic field. If we were 
solving for the magnetic field instead, a similar effect 
would be seen. The magnetic field would rise in the disk, 
but is kept frozen at the boundary, thus building up a 



high magnetic pressure with equally damaging effects for 
the simulation. 

By avoiding the inner boundary altogether, such an 
effect does not occur. However, other problems arise as 
the advection on the very inner disk (s<0.4) happens on 
tight circular trajectories that are poorly resolved in the 
Cartesian geometry near the center of the grid. As they 
behave like highly localized vortices, in most cases the 
explicit shock dissipation terms ensure numerical stabil- 
ity. But for cases with stronger turbulence (model D, Sec 
4.4), a model without an inner boundary could not be 
treated. 

As we did for the models with an inner boundary, 
we check the dependency of the turbulent stresses with 
sound speed by computing versions of models B and C 
without an inner boundary. As for model A2, these mod- 
els B2 and C2 behave quite similarly to models B and C, 
saturating at roughly the same stresses. 

4.4. Radially varying field - models D and E 

With the constant vertical field, it is seen that the inner 
disk does not go turbulent. This is due to the fact that in 
this rapidly advecting region, the growth of the MRI is 
numerically damped. 

In view of this, we also compute models with a radi- 
ally varying z-field 

B = ^^Vw*, (20) 

such that j Balbus-Hawley wavelengths (Eq. fT7\ ) are re- 
solved through the vertical extent L z of the box at any 
radius. We computed a model with j = 4 (model D), 
and another version, with a weaker field (model Dw) 
with j = 16, so that unstable wavelengths are reason- 
ably well resolved. These fields can make the whole disk 
go turbulent, but they grow so strong in the inner disk 
of model D that a more pronounced temperature profile 
(cj T =2) had to be used to avoid the magnetic field from go- 
ing superthermal. For model D, such a setup has plasma 
(6=20 at s; nt and 120 at s ext . Model Dw has a field with a 
strength more similar to model A, and therefore would 
not need this fix. However, to allow a comparison with 
model D, we also used this steeper temperature gradient, 
thus having a plasma /5= 300 at s- mt and 1850 at s ex t ■ 

As this setup has an inverse /5 profile compared to 
the models with a constant vertical field, the turbulence 
starts from the inner disk instead of the outer. Also, as the 
magnetic field is stronger, wavelengths of faster growth 
rate are better resolved and the turbulence saturates after 
a few orbits. Due to the strong stresses, we had to raise 
shock resistivity to 2 instead of 1 as used before. 

The stresses also saturate at higher values, 0.08 for the 
Maxwell stress, and 0.02 for the Reynolds stress (Fig. flT) 
for model D. The weaker field yields a total alpha vis- 
cosity around 2 x 10~ 2 (a M =0.017, a R =0.005, see Fig. Ink 
and Table 2). These values are at least one order of magni- 
tude higher than the ones obtained with a constant field. 
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Fig. 7. Time evolution of the turbulence for different sound speed profiles. The strength of the angular momentum 
transport differs with sound speed. The increase in stress observed when the sound speed is raised to c Sg =0.20 is 
dramatic. Time is quoted in orbits at sq. 



The radial structure of the alpha value, plotted in Fig. [T2j 
reveals that the stresses follow the radial profile of 1 //3, 
being stronger in the inner disk. The same was seen in 
the simulations with a constant field, where the stresses 
were stronger in the more magnetic outer disk. 

The stresses in this model are so high that the alpha 
viscosity in the saturated state is always of the order 
10 _1 , reaching w0.5 in the more magnetized inner disk 
(Fig. [12} ■ That means that the inner disk approaches mag- 
netically dominated values for plasma /3 as the turbu- 
lence saturates. The global average of plasma /3 in the sat- 
urated state is 4, although it did not reach superthermal 
values (<1) during the course of the simulation. Model 
Dw is milder, having a total alpha value of the order 10 -2 , 
reaching a maximum of ~0.08 in the more magnetized 
inner disk. 

The high stresses compared to the cases with constant 
vertical field stem from the initially stronger magnetic 
field, not from a pressure effect. Indeed, in this setup, the 
temperature profile is steeper, but it never rises above 
the sound speed of model C, with c s „ = 0.2 and power 
law of exponent 0.5. However, we still expect the pres- 
sure effect seen on models ABC to be present. To check 
the behavior of this setup with the imposed temperature 
profile, we compute another model (model E) with the 
same initial condition for the magnetic field as model D, 



but with a hotter temperature, with c Sg = 0.2. This setup 
has a plasma beta value of 80 in the inner disk and 450 
in the outer. A weaker version with plasma beta value 
of 1200 in the inner disk and 7400 (model Ew) was also 
computed, for comparison with model Dw. The results 
(Fig [Tib) are similar to models ABC: The stresses are in- 
deed larger than those of model D and Dw, but the al- 
pha viscosity values are smaller. Also as seen on models 
ABC, the kinetic alpha did not change appreciably, but 
the magnetic alpha was reduced by approximately a fac- 
tor of 2. 

The ratio of stresses varied as well, as seen in models 
ABC. It is around 4 for model D, but 2 or less for model 
E. It is around 3 for model Dw, but less than 2 for model 
Ew. 

A final note on the temperature effect; we see that, 
like alpha, the turbulent plasma beta parameter j6f does 
not scale with pressure. As the rise in the turbulent mag- 
netic energy is outpaced by the growth of the pressure, 
jit increases with increasing temperature. 

4.5. Constant azimuthal field - models F and G 

Analytical treatment (Balbus & Hawley 1992, Ogilvie & 
Pringle 1996) and numerical simulations (Hawley 2000, 
Papaloizou & Nelson 2003) show a wealth of evidence 
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Fig. 8. Evolution of the alpha parameters for different sound speed profiles. This quantity, that measures the strength 
of viscosity through the parametrization Vt = ac s H, decreases with sound speed for the magnetic stresses but stays 
constant for hydrodynamic stresses. As seen before, the stresses actually increase with increasing temperature, so this 
decrease of the alpha parameter is due to the stresses increasing less rapidly than the pressure. Time is quoted in 
orbits at sq. 



that the MRI also exists for a purely toroidal field. In this 
case, waves of the form exp[i(m<p — cot)], where m=k^s 
is the azimuthal wavenumber, are excited. The maxi- 
mum growth rates are similar to those observed in purely 
vertical fields, but reached at much smaller azimuthal 
wavenumbers. 

For an azimuthal field, the maximum growth rate oc- 
curs at the wavenumber (Balbus & Hawley 1998) 



m r 



15 ns 



(21) 



Such wavelengths are now resolved in the xy plane, in- 
stead of in the vertical direction. Without the severe con- 
straint of fitting unstable wavelengths in the tiny ver- 
tical scale height of the disk, the azimuthal field can 
be set at much stronger values than those used in the 
vertical cases. The only constraint is that we keep the 
field subthermal. With a temperature gradient of q T =l 
and c So =0.05, a constant azimuthal field of Bq=3 x 10~ 2 
(model F) corresponds to plasma /5 of 12 at s=0.4, 5.5 at 
So and 2 at s=2.5. The wavenumbers are m max =65 at So, 
102 at s=0.4 and 41 at the outer boundary. A hotter ver- 
sion, with c So =0.20, yielding plasma /3=220 at s=0.4, 90 at 
So and 35 at s=2.5 (model G), was also computed for com- 
parison. 

Following the time evolution of model F, we see that 
the turbulence actually saturates at £=200 at s=2.0 (~ 11 
orbits), f=300 at s=1.5 (« 30 orbits) and is still growing 



linearly at s=1.0 at the end of the simulation at t=628 
(100 orbits). The global average (fig.H3l, and a space-time 
(s, i) inspection (fig. fl4)l of the stresses reveals that af- 
ter reaching saturation at £=250 (40 orbits at Sq), a steady 
state is maintained for 20 orbits, after which the turbu- 
lence starts decaying slowly. But as a small growth is ob- 
served near the end of the simulation, it is not clear if this 
decay would continue to zero or if it constitutes just a un- 
usually long fluctuation of the turbulence. Moreover, the 
magnetic and kinetic energies (fig. [131 show no signs of 
decaying. 

On the right panel of fig. [TU we plot the total alpha 
viscosity parameter a M + oc R . Curiously, it does not show 
the decaying effect seen on the stresses, implying that 
a decrease in gas pressure accompanied the decrease in 
Maxwell stress. Such a behavior is expected, since a neg- 
ative density gradient is arising from the accretion pro- 
cess. Therefore, a depleted outer disk has a larger value 
of alpha viscosity for the same Maxwell stress. 

The global average yields R s ^=(0.7±0.1)xl0" 5 , 
M s ^=(2.9±0.4)xl0" 5 and total alpha viscosity 
a=(1.3±0.1)xl0~ 2 . 

Model G shows Maxwell stresses that go further to- 
wards the inner disk, and Reynolds stresses that extend 
to inside s=0.4 as seen in Fig. [TU The outer disk attains 
a turbulent state first, at £=100 (15 orbits), and the inner 
disk at £=200 (30 orbits). The turbulent alpha parameter 
peaks at 8 x 10" 3 (a M =6 x 10~ 3 , a R =2 x 10~ 3 ) at 30 orbits 





Fig. 9. Density contours at selected planes x=0, y=0, and z=0 on saturated turbulent state for sound speed profiles of 
c S g = 0.05 (model A, left panel) and c Sg — 0.20 (model C, right panel). The color code is the same for both figures. The 
stronger stresses for the hotter case lead to a much more effective turbulent viscosity, as seen from the steep density 
profile resulting from accretion. The vertical planes are stretched to show more detail than the correct aspect ratio 
would allow. Movies of these simulations can be found at http : //www. astro .uu. se/~wlyra/planet .html 



and starts a long decay until leveling after other 30 orbits 
at ftM = 3 x 10~ 3 and a^=7.5 x 10~ 4 . The global averages 
are R S( P=(5±2)xlO- 5 l M s< f=(1.8±0.6)xl0" 4 and total al- 
pha viscosity a=(4.2±1.5) x 10~ 3 . 

The referee, Dr. Ulf Torkelsson, pointed out to us that 
the non-zero tension of the constant azimuthal field, 

}T l {B- V)B = -^ _1 B§s _1 S, 

leads to an increase in the centripetal force that is not 
taken into account by Eq. (12) , so the models with az- 
imuthal flux are not started in strict magnetohydrostat- 
ical equilibrium. We therefore performed a set of 2D tests 
to assess how this out-of-equilibrium initial condition 
could modify our results. First we ran a 2D disk without 
noise in the velocity field, so although a magnetic field is 
present, the spectrum of wavelengths is not excited. The 
departure from equilibrium launches a sound wave start- 
ing from the inner disk and propagating outwards. At 
time f =30 (~5 orbits) in model F the sound wave reaches 
the outer boundary and is damped by the buffer zone. 
After that, the oscillations slowly damp through the next 
orbits as the disk settles into centrifugal equilibrium be- 
tween gravity, thermal pressure, and magnetic tension 
forces. We followed the evolution until time £=90. At this 
time, the amplitude of the perturbation dropped to 3% of 
the initial density and 1% of the reference sound speed 

CsO- 

In model G, with a sound speed 4 times faster, the 
sound wave reaches the outer boundary much earlier, 
at time £=12 (about 2 orbits). By including noise we see 
the same results, so we conclude that in a 2D case, even 
though the non-vanishing magnetic tension leads to an 
out of equilibrium initial configuration, the discrepancy 



is slight and the system quickly relaxes in a timescale that 
is much smaller than the time the MRI takes to saturate 
(20 orbits). 

5. Disks with solid boulders 

Having presented the gaseous disk models, we now pro- 
ceed to study the behavior of solid boulders inserted in 
these disks we have constructed. Meter-sized boulders 
are an important step towards kilometer-sized planetes- 
imals. They are also interesting from a gas-dynamical 
point of view because they are only marginally coupled 
to the gas (on approximately a Keplerian shear time- 
scale) and can thus experience concentrations in vor- 
tices and transient gas high pressures (Barge & Sommeria 
1995, Fromang & Nelson 2005, Johansen et al. 2006) 

Our models are usually evolved for fa 75 orbits at So 
before we add the particles, to allow for the turbulence to 
develop and saturate. A large number of particles (10 6 ) 
is used, which allows us to trace the swarm of parti- 
cles onto the grid as a density field. The initial condition 
is such that the particles are concentrated in an annu- 
lus of constant bulk density p p , ranging from Si n t to s ex t , 
with a solids-to-gas ratio of 0.01, typical of the interstellar 
medium. Their velocity is initially the Keplerian angular 
velocity for their radial location. 

As boundary conditions, we do not allow the particles 
to leave the computational domain as they drift inwards 
due to gas drag from the slightly sub-Keplerian gas. 
Instead, a particle that crosses the inner radius Si nt will 
be relocated to the outer radius s ex t , where it will reap- 
pear at the same azimuthal location, thus mimicking pe- 
riodic boundary conditions. Its velocity, however, will be 
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Fig. 10. Radial density profiles for models A and C. The lower panels show space-time plots of the vertically and 
azimuthally averaged density Model A shows what seems to be a variability around a constant value, model C has 
developed a smoother radial density gradient. The lower panel shows the two density profiles at the end of the 
simulation. 



reset to the Keplerian angular velocity at s ex t . A particle 
that tries to cross the outer radius will simply have its 
velocity reset to Keplerian without changing its radial or 
azimuthal location. 

We include particles at the later stages of models A 
and D. For model A, where the inner disk does not go 
turbulent, it was noticed that by allowing the particles 
to move through all the radial range, they eventually got 
trapped in the several local density maxima between the 
concentric rings that the inner disk breaks into. As such 
a loss of particles is undesirable, we keep them where 
the turbulence is saturated by setting the inner radius 
of the boundary conditions for particles outwards of sq. 
In model D, where the Balbus-Hawley wavelength is re- 
solved at all radii, such correction is not needed. We dis- 
cuss the results of model D first. 



5.1. Particles on model D 

The particles soon fall to the disk midplane due to the 
vertical gravity. At the same time they get trapped in high 
pressure regions in the turbulent flow due to the drag 
force (Klahr & Lin 2001, Johansen et al. 2006). In Fig. [ll 
we plot a slice of the bulk density of solids profile 10 or- 
bits after the insertion of the particles into the simulation. 
A pile-up of solids in the inner disk is seen to have oc- 
curred, because particles have concentrated in a pressure 
maximum in the gas (we discuss this further in Sect. 15. 2\ 
see also Fig. ll8|l . 

As discussed by Johansen & Klahr (2005), while solid 
particles are pulled towards the midplane by the stellar 
gravity, turbulent motions stir them up again. A sedi- 
mentary layer in equilibrium between turbulent diffu- 
sion and gravitational settling is formed. The thickness 



16 



10" 



1(T 



Lyra et al.: Disk-in-a-box 

Stresses - Model D Stresses - Model E 

1(T 



10" 




10 20 30 40 50 60 

//(27i^ 1 ) 

Stresses - Model Dw 



10" 



10" 



10- 




10" 



10- 



20 40 60 80 100 

t/(2ltQ.g) 



10" 



10" 




10 20 30 40 50 

?/(2ji£2„ 1 ) 

Stresses - Model Ew 



8 



10" 



10- 



10- 




20 40 60 80 100 

f/(27ifl„ 1 ) 



Fig. 11. Time evolution of the turbulent stresses for models D and E (upper panels) and Dw and Ew (lower panels). 
As compared with model A, the stresses saturate at a much earlier time. The ratio of Maxwell to Reynolds stress is 
around three for model D, and around two for the hotter model E, approaching one at the end of the simulation. 
Notice that as seen in models ABC, the stresses are bigger for the hotter model, although the alpha viscosity value is 
smaller. Time is quoted in orbits at Sq . 



of this layer is therefore a measurement of the turbulent 
diffusion acting on the solid particles. 

Under the influence of gravity, the solids settle with 
a profile similar to the one generated by a pressure force 
(Dubrulle et al. 1995) 

z 2 

lnp p (s,z) =lnp p (s,z = 0)-^2. (22) 

By comparing this profile with the analytical expres- 
sion for a pressureless fluid under diffusion, gas drag and 
vertical gravity (Johansen & Klahr, 2005) 

mpp = lnp p (s,z = 0) {- \ g z dz, (23) 

d¥> J 

and recalling that g z = — H 2 z, we have 

D z t] = n 2 H 2 T/ (24) 

From Eq. ll22l , we see that the scale height of the solids is 
the vertical distance in which the bulk density falls by a 
factor 1/ \fe ps 0.6 relative to the value at midplane. We 



plot in Fig. ITBb the bulk density normalized by its value 
in the midplane. In this figure, the quantity plotted is in 
fact identical to the exponential term in Eq. \22\ . Where 
it reaches 0.6, the vertical distance z gives the diffusion 
scale height Hp. 



We fit the points where the exponential term equals 
0.6 with a power law Hp = ar" . A linear regression in 
logarithm yields a = 0.042 and n = 0.97, with an rms of 
0.04. This translates into a diffusion coefficient (Eq. l24l ) 
of D^' pa 1.7 x 10 -3 s -1 . As this model has a sound speed 
profile c s =0.1 s _1 , the diffusion coefficient in dimension- 
less units corresponds to <?W = D^c^ 2 Cl = 0.17s -15 , or 
0.14 if globally averaged. The rms of 0.04 in the logarithm 
fit yields an uncertainty of 0.01 in this global average. As 
the total alpha viscosity is 0.112±0.003, the globally aver- 
aged vertical Schmidt number, i.e., the strength of viscos- 
ity when compared to vertical diffusion, is 0.78±0.06. 




Fig. 12. Same as Fig.[6j but for models D (left panel) and Dw (right panel). The times corresponding to the snapshots 
are indicated in the legends, as well as the time average. The inner disk is considerably more turbulent than the outer 
parts. 
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Fig. 13. Same as Fig. [5] (constant vertical field), but for model F (constant azimuthal field). 



5.2. Particles in model A 

For model A, the alpha viscosity (a ~ 10~ 3 ) is much 
lower than in model D (a ~ 10 _1 ), so according to 
Eq. l|24t and assuming that the diffusion coefficient is of 
the same order of the turbulent viscosity, we expect the 
sedimentary layer of solid particles to have a scale height 
of Hp = u/€L 2 Tf « 0.03H. At s ext , the gas scale height H 



equals 0.08 for c So = 0.05, then H p = 2.5 x 10~ 3 . With a 
grid resolution Az = 0.02, this layer will not be resolved. 
It means that the interpolation of the particle density back 
to the grid will not allow for a grid-based measurement 
of the diffusion acting in the turbulent layer as we did for 
model D. 

But as the particles are Lagrangian, we can plot their 
real positions and trace the scale height of solids in a 
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Fig. 14. Space-time diagram of the turbulent stresses and alpha viscosity for model F (upper panels) and model G 
(lower panels). Time is quoted in orbits at So=1.0. Saturation is reached at 40 orbits, but after 60 orbits the stresses 
seem to start a slow decay. The alpha viscosity appears constant, due to a similar decay in gas pressure in the outer 
disk, that starts to deplete as the resulting accretion builds a negative density gradient. Model G behaves similarly, 
but with higher stresses, lower alpha viscosity and being turbulent further inside. 



grid-independent way. The diffusion process operates in 
much the same way, nearly independent of grid resolu- 
tion, since the large scale velocities of the gas are well 
resolved. In order to do this, we define 128 bins in the 
radial direction and measure the individual vertical po- 
sitions of the swarm of particles with respect to the mid- 
plane of the disk within these bins. The standard devia- 
tion z rms of the vertical positions of particles with respect 
to the midplane in each bin immediately gives the scale 
height of the sedimentary layer. The result of this process 
is shown in Fig. \l6k , where we average z rms as measured 
on 17 snapshots, from orbits 4 to 20 at so after the inser- 
tion of the particles. The initial time is chosen at 4 orbits 
because it is the time it takes for the drag force to cou- 
ple the particles to the gas at s ext , thus making sure that 
the sedimentary layer is in equilibrium between gravita- 
tional settling and turbulent diffusion through the whole 
radial extent of the disk. 

The dashed line in Fig. [I6h represents the power 
law fit Hp=ar n to the measured scale height. It yields 
fl=0.003 and n=2.48. The rms of the logarithmic fit is 



0.09. In Fig. [16b we plot the resulting diffusion coefficient 
Dz =Cl 2 HpTf. It can be approximated by a power law of 
« 9 x 10~ 6 s 2 . In dimensionless units it corresponds to 
£( f )=3.6 x 10 _3 s -1 - 5 , or 0.007 if globally averaged. The 
uncertainty is 0.001. This behavior of the turbulent diffu- 
sion that acts on solids is quite similar to the one shown 
by the turbulent viscosity arising from the stresses on the 
gas phase (dot-dashed line in Fig. [16b). Checking the ra- 
dial dependency of the Schmidt number (Fig.H6b), we see 
that it is of the order of unity all over the radial domain. 
A slight trend is seen towards smaller Schmidt numbers 
in the outer disk, but it never gets below 0.6. 

In Fig. [17] we explore this radial dependency in more 
detail. The figure shows the vertical and azimuthal aver- 
age of the turbulent viscosity Vt, turbulent diffusion Df 
and their ratio Sc, the Schmidt number, as a function of 
radial position s. On the time axis we show time in or- 
bits at So since the beginning of the simulations. At 71 
orbits the particles are inserted, and quickly fall to the 
midplane. As seen from the diffusion map, the particles 
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Fig. 15. Vertical slice of the bulk density of solid particles (upper) and the same quantity normalized by the midplane 
density. A midplane layer forms in equilibrium between sedimentation and turbulent diffusion. The scale height of 
this layer follows a linear dependence with radius. 



in the innermost radii quickly sediment, settling in a dif- 
fusive equilibrium in less than two orbits. After four or- 
bits, t = 75, the last radius achieves diffusion equilibrium 
and the situation becomes statistically unchanged until 
the end of the simulation. 

The gas-phase viscosity and the solid-phase diffusion 
are similar in average, but the diffusion is seen to fluc- 
tuate more. Transient patches of high diffusivity are seen 
to live for some orbits at a constant radial location be- 
fore decaying. The turbulent viscosity, in turn, appears 
much smoother in time. The resulting Schmidt number is 
shown in greyscale in the right panel of Fig. [17| Due to 
the high variability of diffusion, some short lived bright 
areas of Sc>5 are seen, but overall the Schmidt number 
is around 1 throughout the space-time domain. 

In Fig. [18] we show contours of the gas density (left 
panel) and the bulk density of solids (right panel) in a 
snapshot taken after 20 orbits after the insertion of the 
particles. A correlation is seen as the solids show large 
concentration at areas of high gas density. Initially, the 
density increases linearly as the particles sediment to- 
wards the midplane. As seen before, the sedimentation 
is complete at the outer radius after ~4 orbits. After that, 
the growth is only due to the particles being concentrated 
in transient gas high pressures. 

The average bulk density is quite low (p p =0.003, or 
6.0 x 10 _11 kgm -3 in physical units, see Table 1), and 
several areas devoid of particles are seen in the disk. 
However, the overdensities observed as bright clumps in 
the snapshot are several standard deviations above av- 
erage. By plotting the maximum solid density (which is 



roughly the solids-to-gas ratio) throughout the simula- 
tion (Fig. Il9fr , we see that its value is usually around 30, 
but it can reach values as high as 85. Such a behavior was 
seen in shearing box simulations by Johansen et al. (2006), 
who report local enhancements of the solids-to-gas ratio 
by a factor of 100, also pointing that such concentrations 
are gravitationally unstable, thus being able to collapse 
to form km-sized bodies. 

As a control, we also simulated the settling for a non- 
turbulent, purely laminar, unperturbed disk. Without 
high pressure regions, the particles simply sedimented 
towards the midplane forming a thin homogeneous layer 
of solids. 

6. Summary and conclusions 

We have considered MHD models of global Keplerian 
disks in Cartesian grids. These disk-in-a-box models are 
able to develop and sustain MHD turbulence, in good 
agreement with published results achieved with cylindri- 
cal codes and shearing boxes. In this first article of the 
series, we investigated the dependence of the MRI with 
disk scale height and the dynamics of solid boulders in 
the global hydromagnetic turbulence. 

As a numerical solver we have used the PENCIL 
CODE. This finite-difference code solves the non- 
conservative form of the dynamical equations using 
sixth-order spatial derivatives, achieving spatial resolu- 
tion that approaches that of spectral methods. The nu- 
merical scheme is stabilized by using hyperdissipation 
and shock dissipation terms, which enter as free parame- 




Fig. 16. Time-averages of the dispersion of vertical positions of particles w.r.t. the disk midplane, revealing the scale 
height of the sedimentary layer that forms due to the equilibrium between sedimentation and turbulent diffusion 
for model A. The dotted line is a power law fit, yielding an exponent of « 2.5. This radial profile is a time average 
between orbits 4 and 20 (see text). The resulting diffusion coefficient resulting from this scale height is shown as solid 
line in the middle panel. The time averaged gas-phase viscosity is shown in dot-dashed line for comparison. It is seen 
that these two coefficients have similar strength. The right panel shows the Schmidt number of the flow. It is indeed 
close to unity through the radial domain. 
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Fig. 17. Azimuthally and vertically averaged turbulent viscosity Vt, turbulent diffusion Dt and the resulting Schmidt 
number Sc as a function of radial position and time for model A. Time is quoted in orbits at So- Viscosity and diffusion 
are shown in the same units and color-code. Some localized overdensities in diffusion last for some orbits, while 
viscosity shows a smoother evolution in time. The quantities have approximately the same strength, as the Schmidt 
number is overall around unity, and seldom greater than 5. 



ters in the dynamical equations. The effect of hyperdissi- 
pation is to quench unstable modes in the small scales of 
the grid, while affecting the large scale motion as little as 
possible, whereas shock dissipation is invoked to smear 
out large divergences in the flow field. We choose these 
parameters by performing series of 2D gap opening sim- 
ulations with a Jupiter mass planet and comparing them 
with a higher resolution calculation without explicit dis- 
sipative terms. 

We find evidence that the turbulence generated by 
the magnetorotational instability grows with the thermal 
pressure. The turbulent stresses depend on thermal pres- 
sure obeying a power law of 0.24 ± 0.03, compatible with 
the value of 0.25 found in shearing box calculations by 
Sano et al. (2004). We extend this result to a global disk 
showing that the rise in pressure increases the turbulent 
stresses, thus raising the angular momentum transport 



(and therefore the mass accretion rate) although the al- 
pha viscosity value drops. 

We also notice two curious effects. First, the domi- 
nance of the radial component of the turbulent kinetic 
energy increases with temperature. The percentage of the 
total kinetic energy stored in the radial component is 40% 
for the cold model A, and 60% for the hotter model C. 
Second, the ratio of stresses —M S( t'/R s 'l' diminished with 
increasing temperature. It is 5 for model A, and just 1.3 
for model C. The same is seen in the model without inner 
boundary, where the ratio is 6.5 for model A2 and very 
close to 1 for model C2. This effect is unexpected since it 
is believed that the shear parameter alone controls the ra- 
tio of stresses (Pessah et al. 2006, Ogilvie & Pringle 1996). 
From the shearing box data of Sano et al. (2004) the stress 
ratio seems to be constant with temperature. 
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Fig. 18. Density contours at midplane of the gas and solid phases of the disk (left and right panel, respectively). The 
snapshots were taken at 20 orbits at sq after the insertion of the particles. The color code for the solid phase is selected 
to represent 2 sigma (0.25) above the average bulk density of 0.03. The bright areas are saturated as the maximum den- 
sity reaches as far as 85. A movie of this simulation can be found at http : //www . astro .uu.se/ ~wlyra/planet . html 
A correlation with gas density is seen, since the bright clumps of solids correspond to pressure maxima, i.e, areas of 
high gas density. 



One explanation could be that, according to Eq. [T2j 
the angular velocity is sub-Keplerian and the increasing 
effects of pressure from the colder to the hotter mod- 
els modifies the shear. Quantitatively, however, one sees 
that the pressure correction is too small to account for 
the decrease in the stress ratio and, more importantly, 
would have the opposite effect. According to the lin- 
earized equations for the evolution of the turbulent fluc- 
tuations, the Maxwell stress couples with shear qCl, and 
the Reynolds stress couples with the large scale vortic- 



ity w = (2 - q)Cl (Balbus & Hawley 1998), where q = 
— 31nQ/31nr is the shear rate. The pressure-corrected 
angular velocity of the gas can be approximated from Eq. 
CE} as 

n~n K (i-//), (25) 

where rj = (l/2)(31nP/31nr)(H/s) 2 > is a parame- 
ter often used to parameterize the strength of the global 
pressure gradient (see e.g. Nakagawa et al. 1986). Typical 
values of q lie between 0.001 and 0.1. 
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Fig. 19. Maximum bulk density of solids, in units of the 
mean gas density as a function of time for model A. Time 
is quoted in orbits at Sq. The maximum density rises as 
the particles sediment towards the midplane. After the 
sedimentation that lasts for four orbits, the particles are 
coupled with the gas and trapping in transient gas high 
pressures raise their maximum density well above aver- 
age. The maximum density is usually around 30, but be- 
tween orbits 12 and 13, it reached values as high as 85. 



The reduction of both the angular frequency and 
shear rate should reduce the Maxwell stress. Our sim- 
ulations show the opposite, with the Maxwell stress in- 
creasing as the pressure is raised. Regarding the stress 
ratio, reducing the shear increases this quantity since 
the Reynolds stress falls faster than the Maxwell stress 
due to the stabilizing effect of the growing vorticity 
(Abramowicz et al., 1996). Once again, we see the oppo- 
site effect. 

As most of the analysis of turbulent thin accretion 
disks have focused on locally isothermal simulations us- 
ing c s « 0.05, changing the field configuration while 
keeping the temperature constant, such behavior has 
been largely overlooked. Although the disk temperatures 
considered in this case are quite extreme for disks around 
T-Tauri stars, circumplanetary disks are thought to be 
rather thick (Klahr & Kley 2006) and therefore the evo- 
lution of the MRI in such disks is expected to be more 
similar to the hotter cases considered in this paper (mod- 
els CEG) than the colder ones. 

We investigated the effect of an inner boundary in 
the evolution and outcome of the turbulence. By using a 
Cartesian grid, an inner boundary can be discarded pro- 
vided we smooth the gravitational potential to avoid a 
singularity in the flow. Models without an inner bound- 
ary do not show the spurious build-up of magnetic pres- 
sure and Reynolds stress seen in the models with bound- 
aries, while the global stresses and alpha viscosities are 
similar in the two cases. 

In treating the solids, we make use of a large num- 
ber of particles, which allows us to effectively map the 
particles back into the grid as a density field without 



using fluid approaches. We monitor the settling of the 
particles toward the midplane and the formation of a 
sedimentary layer when the solids are subject to gas 
drag and the gravity from the central object. The effec- 
tive diffusion provided by the turbulence prevents fur- 
ther settling of solids, in accordance with the results of 
Johansen & Klahr (2005). By having the global disk per- 
spective, we could measure the radial dependence of the 
diffusion scale height of the solid component. The mea- 
sured scale heights imply turbulent vertical diffusion co- 
efficients with globally averaged Schmidt numbers of 
1.0±0.2 for model A (a « 10~ 3 ) and 0.78±0.06 for model 
D(a « lO" 1 ). 

We conclude that the models presented in this first 
paper of the series are capable of sustaining turbulence 
and are adequately suited for further studies of planet 
formation. Future papers will present studies of thermo- 
dynamics and radiative transfer in the evolution of the 
turbulent stresses, planet-planet and planet-disk interac- 
tion, the effect of stratification and the dynamics of dead 
zones. 
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Appendix A: Anisotropic hyperdissipation 

Hyperdissipation is used to quench unstable modes at 
the grid scale, therefore being intrinsically resolution- 
dependent. Because of this, isotropic dissipation only 
gives equal dissipation in all spatial directions if Ax = 
Ay = Az, i.e., if the cells are cubic. For non-cubic cells, 
anisotropic dissipation is required as different directions 
may be better /worse sampled, thus needing less /more 
numerical smoothing. Such a generalization is straight- 
forward. We notice that hyper diffusion works as a con- 
servative term in the continuity equation such that 

f D (p) = V-J, (A.l) 

where 3 = D$'V 5 p is the mass flux due to hyperdif- 
fusion. For simplicity, we will drop the subscripts "3" 
from the coefficients hereafter. This formulation reduces 
to the usual sixth-order hyperdiffusion under the condi- 
tion that D is constant. Generalizing it to three dimen- 
sions simply involves replacing this mass flow by 
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,D Z 
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so that different diffusion operates in different directions. 
Since D x , Dy and D z are constants, the divergent of this 
vector is 



V-J = D x 



L dx b 



D 
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The formulation for resistivity is strictly the same. For 
viscosity it also assumes the same form if we consider 

a simple nth order rate of strain tensor operator sjj- = 
Appendix B: Shocks 

Shock viscosity is taken to be proportional to posi- 
tive flow convergence, maximum over three zones, and 
smoothed to second order, 

Cv = v s h /max[(-V ■ »)+]\ [min(Ax, Ay, Az)] 2 , (B.l) 

where v s ^ is a constant defining the strength of the shock 
viscosity, usually around unity. We refer to it as the shock 
viscosity coefficient (Haugen et al. 2004). In the equation 
of motion it takes the form of a bulk viscosity so that now 
the stress tensor contains 

Tij= [...]+ pC v SijV-u, (B.2) 

where [...] refers to the (hyper) viscous terms described 
in Appendix[A] The acceleration due to shock viscosity is 
therefore 

fv(u,p) = |0 _1 V ■ T 

= ^[V(V- M ) + (Vlnp + Vln^,) V ■ w] . 

(B.3) 

Such a viscosity scheme ensures that energy is dissi- 
pated in regions of the flow where shocks occur, whereas 
more quiescent regions are left untouched. The formula- 
tions for shock diffusion and shock resistivity are similar, 
yielding 

/ D (p) = £ D (v 2 p + Vln£ D - Vp), (B.4) 

and 

/, (A) = £, ( V 2 A + V In V ■ a) , (B.5) 

where £rj and are analogous to £ v in Eq. l|B.lt , contain- 
ing their respective shock diffusion and shock resistivity 
coefficients D sh and w sh . 
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